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Abstract 

The modeling of EEG with time-varying autoregressive (TVAR) models in which the 
parameter evolution is constrained to be a linear combination of basis functions is discussed. 
As a suggestion how to overcome the difficulties inherent in using generic basis functions and 
separately adjustable coefficients, a method is presented in which the basis is constructed 
so that the expected parameter evolutions can be optimally approximated by the basis. The 
method is then applied to the estimation of event related synchronization changes in the EEG. 
The results show that a properly selected basis enables effective estimation of the parameter 
evolution if the process to be estimated is compatible with the learning set with which the 
basis were constructed. 

1 Introduction 

Autoregressive (AR) time series models have been successfully applied to the analysis of stationary 
EEG epochs since the end of the 1960's [6], [29], [7]. The ends to the modeling have been diverse, 
such as simulation [28], spectral analysis [7], classification of sleep stages [11], estimation of the 
depth of anaesthesia [4], data compression and transient detection [24], prediction of epileptic 
seizures [27] and estimation of information flow between different parts of the brain [8]. The 
usual way of dealing with the modeling of nonstationary EEG is to use either fixed or adaptive 
segmentation [3], [16]. The EEG is then treated as a stationary process within these segments. This 
approach has been proven to work satisfactorily in many cases, such as in sleep when the EEG 
typically changes from one state to another. Another such a case is when the EEG exhibits slow 
trends. An example is the slowly changing effects of medication. 

In some cases the rate of change of EEG is such that if we use segmentation the segments 
turn out to be so short that the estimates are useless due to their small sample properties, see 
e.g. [25] for approximations for the confidence intervals of AR parameters and the corresponding 
spectrum estimates. One obvious solution is to allow the AR parameters to be time varying. Such 
models are called TVAR models. However, we can not assume that the parameters are arbitrarily 
time- varying since the parameter estimation problem would be severely underdetermined with an 
infinity of parameter evolutions that yield vanishing prediction errors. Such parameter evolutions 
are clearly meaningless. One approach to make the estimates meaningful is to assume that the 
parameter evolution is either slow, smooth or both. 

There are two main classes of methods to fulfill these constraints. The first is to use recursive 
estimation of the time- varying predictor evolution and the second is to constrain the parameter 
evolutions to be linear combinations of some basis functions with appropriate properties. These 
have been called the stochastic and deterministic regression approaches, respectively [9]. We 
shall concentrate here on the deterministic regression approach that we shall call the least squares 
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TVAR (TVARLS) scheme, since the actual estimates are formed so that the norm of the residuals 
is minimized. 

The TVARLS scheme was first introduced in 1970 [26]. The approach was to approximate the 
parameter evolution by truncated Taylor expansion of order M and to solve the weighted least 
squares estimate of the expansion coefficients. This is equivalent to forcing the parameters to be 
polynomials of order M. The idea has been further developed, e.g. in [19], [13], [14], [10] and 
several bases have been proposed. 

In the usual TVARLS scheme the coefficients of the basis for each parameter evolution are 
separately adjustable. If we assume that the process is stable at all times, we can easily run into 
trouble with processes that are temporarily of the narrow band type. The roots of the characteristic 
polynomial of the model will then easily have modulus greater than unity temporarily, which is 
due to the tendency of the parameters to overshoot in the parameter space. Naturally, unstable 
models are not always a serious problem e.g. if the end to the modeling is prediction since the 
predictors corresponding to finite order AR models are always stable. However, if the end is e.g to 
obtain an evolutionary spectrum estimate, this is a major problem. 

In some cases more can be assumed of the process than just smoothness or slow variation. One 
such a case is the event related synchronization/ desynchronization (ERS/ERD) of alpha waves of 
EEG [20], [23]. For example, when a patient has his/her eyes closed, the occipital EEG (primary 
visual cortex) shows high intensity in the 8-12 Hz region (alpha band, synchronization) while the 
opening of the eyes this intensity decreases or even vanishes (desynchronization). In this kind of 
a situation we can assume that the EEG exhibits a more or less rapid transition from a state to 
another and that the transition starts at some time after the visual stimulation. 

While it seems that we had access to an ensemble of realizations by repeating the visual 
stimulation, this is not necessarily the case. It is well known that in many EEG experiments of 
this type the patients are subject to habituation, that is, the responses to stimulations may exhibit 
various kinds of trends. In such cases the treatment of consecutive responses as an ensemble does 
not give us a correct view of the situation since we do not actually have an ensemble. Rather, 
we will obtain some kind of "average" estimates if we e.g. calculated the covariance of these 
"pseudoensembles" . 

In this paper we propose an approach to TVARLS modeling in which prior assumptions on 
the dynamics of the process can be taken into account. This is accomplished by first forming 
a representative set of expectable parameter evolutions, or a learning set. Since the individual 
parameter evolutions in practically all relevant cases are mutually correlated, the learning set is 
constructed so that these correlations are not lost. A low dimensional subspace of the space of 
all expectable parameter evolutions is then extracted in such a way that the evolutions of the 
learning set can be approximated in this subspace with minimal error. This is accomplished via 
the eigendecomposition of the covariance matrix of the learning set. 

The rest of the paper is organized as follows. In Section 2.1 we formulate the conventional 
TVARLS problem. In Section 2.2 we discuss the problems inherent in the use of generic bases 
and the determination of the adjusted bases in general. We also modify prediction equations 
to give a scheme with non-separately determined basis. In Section 3.1 we discuss the formation 
of the learning set for the visual ERD/ERS data. In Section 3.2 we evaluate the performance 
of the method with simulations and in Section 3.3 we apply the method to real ERS data 
and verify that the method is applicable to the ERD/ERS measurements. We also show how 
the evolution estimates could be employed to further estimate e.g. delays in the start of the 
synchronization/desynchronization after the trigger. Finally, in Section 4 we address the problems 
and the potentials of the proposed method and discuss some previous methods that have been 
used in the analysis of ERD/ERS. 
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2 Methods 



2.1 TVAR Model 

The deterministic regression TVAR model for a process (sample) x t can be written as 

v 

x t = ^2a k {t)x t - k + e t , (1) 

k=l 

where e t are the residuals and p is the order of the AR model (predictor). One of the criteria 
for the validity of the model is that the residuals can be modeled as zero mean independent 
random variables with variance p t . Let the observed time interval be t = 1,...,T. In the 
conventional TVARLS scheme the time- varying prediction coefficients (the AR parameters) a k {t) 
are separately constrained to a (usually smooth) subspace S C M. T with dimension M and basis 
<f>i{i)i 02 ) ? • • • j 0m(O where the choice <j)\{t) = 1 is usually made to include the stationary model 
as a special case. The time- varying parameters are thus of the form 

M 

a k(t) = ^ C mk(/>m(t) • (2) 
m=l 

Denote the p-order linear predictor of x t by x t with 

v 



k=l 
M p 

= 1 ^ 1 ^2c mk <j) m (t)x t - k (3) 



m=l k=l 

so that x t = x t + e t . 

Instead of the usual approach to focus on the normal equations we will employ the regression 
formalism since the modification to the TVARLS scheme that we propose in this paper is more 
understandable in the latter formalism. To obtain the nonwindowed prediction equations we write 
(1) in matrix form for t = p + 1, . . . ,T and define c = (en, . . . , Ci p , . . . , cmi, • • • , cm p )', A = 
(x p+ i, . . . , xt) 1 # = (ep+ij • • • j e T)' and the regressor matrix 

(0i(p + • • • 0i (p + l)xi • • • (f) M (p + l)x p • • • 0m(p + \ 

; ; ; ; • (-0 

(pi(T)x T -i - - - (/>i(T)x T - P -i • • • 4>m{T)x T -i • • • (/>m(T)x T - p -i J 

where the prime denotes transposition. The prediction equations take now the form 

X = Hc + E . (5) 

If H is nonsingular, the solution c that minimizes E ! E can be expressed in closed form as c = 
(H 1 'H) -1 H' X , that is, the explicit solved form of the normal equations. The estimated parameter 
evolution is then obtained from (2). In this formalism each parameter evolution is estimated 
separately as a linear combination of the basis functions and we have pM free coefficients in (5). 

Several sets of functions have been used as TVARLS bases. These include at least 
ordinary polynomials, Legendre polynomials, sinusoidal (Fourier) bases, discrete prolate spheroidal 
sequences, Walsh and Haar bases [26], [19], [10], [13], [5],[9],[14],[1]. With the exception of the last 
two items, the bases (and their linear combinations) can be described as smooth. 

2.2 Determination of the optimal basis 

It is clear that the ability of the TVAR scheme to model the time- varying characteristics of the 
process is limited by the ability of the chosen set of basis functions to approximate the optimal 
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predictor evolution [1]. Without any prior information on the evolution of the optimal predictor 
coefficients we are not able to treat the problem of optimality of the basis either. 

If the second order statistics of the process were known, the optimal (in the mean square sense) 
time-varying linear p-order predictor for x t could be calculated, e.g. by using the methods given 
in [2] or [15]. Assume that the optimal predictor exists but we have no access to it and use the 
TVARLS scheme instead. If we now have reason to assume e.g. that the evolution of the optimal 
predictor is constant on a time interval, we should employ such a basis that forces the parameter 
evolution to be constant on this interval. Generally speaking, the subspace S should contain only 
such parameter evolutions that are assumed to be expected. Naturally, in most cases this kind of 
prior assumptions can not be made, but when these can be done, they should somehow be taken 
into account. 

Assume further that we have prior estimates on the probability distribution of the predictor at 
times t = 1, . . . ,T and that these distributions are not separable. If we use separately adjustable 
coefficients in (2) we can not introduce prior information on the densities to the estimates easily. 

For these reasons we consider coefficientwise concatenates of the parameters in the sequel. This 
way we can take into account the mutual correlation between the parameters of the predictor. We 
denote these concatenates by <3> £ R K , K = pT 

$ = Ml), . . . , ai(T), . . . , ap(l), . . . , a v (T)) ! . (6) 

The concatenates <3> are also called evolutions although the indexes to <3> are not interpreted as 
previously. 

Let us assume that we have access to a set of evolutions <£ n , n = 1, . . . , TV that can be assumed to 
represent all expectable evolutions <3>* in the following sense. Let G C R K be the subspace spanned 
by {$ n , n = 1, . . . , TV} and Pq be the orthogonal projector onto G. Of all <l> € G the one nearest 
to <l>* with respect to the Euclidean norm is $ = P®$>* . Assume now that ||<£>* — $|| < Si for some 
appropriately small Si. We could thus use a basis for G as the basis for the parameter evolution. 
However, since the dimension TV of G can be large, we can easily run into stability problems when 
we try to estimate the parameter evolution. For this reason we aim next to determine an M- 
dimensional subspace G(M) C G so that the mean of the errors ||<3> n — Pq(m) || 2 3 n = 1, . . . , TV 
is minimal, where Pg(m) 1S the orthogonal projector onto the subspace G(M). 

Denote by T the non-centered sample covariance matrix of the set {<£ n , n = 1, . . . , TV} 

N 

r = TV -1 ®n& n = TV^M' , (7) 

71 = 1 

where the K x TV matrix <5 = (<£>i, . . . , <&/v)- Let 7& and A&, k = 1, . . . , K, be the eigenvectors and 
eigenvalues, respectively, of T, where Ai > A2 > • • • > \k > 0. 

It is well known that the mean error criterion is minimized when G(M ) is spanned by the M 
eigenvectors 7 m corresponding to the M greatest eigenvalues A m , m = 1,...,M [22], [17]. The 
subspace G(M) is called the principal subspace of dimension M and if the set {^n, n= 1, • • • , N} 
was a true ensemble, we would call this approach the sample principal component analysis. We can 
thus use these eigenvectors as the basis for parameter evolution with the modifications discussed 
in Section 2.3. 

The mean approximation error of the learning set in G is 

K 

m=M+l 

The mean normalized approximation error is J^f/trr, where tr denotes the trace of a matrix. If 
the eigenvalues decay sufficiently fast, we may obtain a small normalized approximation error with 
a small value of M. 

Let the maximum approximation error be \\& n — Pe^^nll < ^2 for all n. We have then by 
triangle inequality that 

H** -Pe(M)$* II <S! + S 2 (9) 
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for all expectable parameter evolutions The key problem in the formation of the learning set 
and in the determination of B is to achieve as low a dimension M as possible for a given total 
approximation error 6\ + <fe. The notion of "optimality" now refers to the optimality of B(M) 
with respect to the learning set, not to the expected evolutions. The determination of the optimal 
subspace 0*(M) with respect to the expected evolutions requires that the co variance of these is 
known, which can very rarely be the case. However, if we can assume that Si is small, B(M) can 
be assumed to be a good approximation to B*(M). 

2.3 The modified prediction equations 

The constraining of the (concatenated) parameter evolution to B(M) is equivalent to replacing 
the functions 0 m (£) with the respective sections of j m (t) for each individual parameter and setting 
Cmk = Cm for all k = 1, . . . ,p. 

Let 7 m = (j m (ty, . . . , J m (tyy, m = 1, . . . , M and t = 1, . . . ,T. The predictor coefficients are 
now written in the form 

M 

a k (t) = J2 Wtit) ■ (10) 

m=l 

Thus, instead of the columns of H in (4) we have 

v 

Vrn{t) = Y J l k rrXt)*t-k (11) 

k=l 

as regressors. The prediction estimate is now of the form 

M V 

m=l k=l 
M 

= J2 C m¥m(t) . (12) 
m=l 

The prediction equations are as in (5) but with c = (ci, . . . , cm)' and 

/ <pi(p+l) ••• 9?m(p+1) 



H = 



(13) 



V MT) ••• <Pm(T) 



The dimensionality of the LS problem has now been reduced from pM to M and the estimated 
parameter evolutions are in B(M) by construction. 

It is recommended that the LS solution of c is not done via the normal equations. Instead, the 
reduction of if to a canonical form via an orthogonal transform will lead to a modification of the 
prediction equations, the solution of which is numerically more stable [12]. 



3 Modeling of ERD/ERS 

3.1 Construction of the learning set 

To test the proposed approach we conducted a visual ERD/ERS test using the international 10-20 
electrode system. This was carried out by opening and closing the eyes with an auditory stimulus 
(beep) on 15 second intervals. The subject was a healthy young female. We chose five samples of 
the measured EEG at the transition from desynchronized state to synchronized state (eyes open 
— »■ eyes closed). The samples were low-pass filtered to allow for the decimation of the sampling 
frequency from 250 Hz to 250/4 Hz. The final length of the samples was 960. 

We estimated the AR parameters with p = 6 for both synchronized and desynchronized states 
for all five samples using segments of length 240 at the tails of the samples. The parameters were 
estimated with the modified covariance method [21]. The samples are shown in Fig. la. The 
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Figure 1: a) Five samples of ERS. Vertical scale in microvolts and time 
in seconds. The trigger is set at time t = 0. b) Amplitude spectra of 
AR(6) models of desynchronized (top) and synchronized (bottom)states 
for the five samples. Vertical scale is arbitrary and horizontal is in 
normalized frequency. 



amplitude spectra corresponding the models are shown in Fig. lb. We form the learning set 
first. We assume that the sets of AR(6) parameters describe the EEG in the synchronized (A^) 
and desynchronized (Af) states. We form the parameter evolutions from every Af to every A^ 
assuming that the transition between the two states is smooth. As the model for transitions we use 
sigmoid function for which we use three different slopes and 10 different delays of transition. As 
the length of the parameter evolutions we use T = 240. The parameter processes Aijke(t) G R pxT 
are formed by 

AM*) = A? + (Af - Af)a ke (f) (14) 
where z, j = 1, . . . , 5, k = 1, . . . , 10, £ = 1, 2, 3 and cr^(i) is the sigmoid function 

<r k e(t) = (1 + d 0 exp(d u (t - d 2k )))- 1 (15) 

where d 0 = 1/9, d u = 0.056, d 12 = 0.08, d 13 = 0.12 and d 2k = 68 + S(k - 1). The parameter 
evolutions are monotonous by construction. The extreme cases of the sigmoids o~ke(t) are shown in 
Fig. 2a. The learning set thus contains N = 5 • 5 • 10 • 3 = 750 vectors <3> n of length K = pT = 1440, 
that are the componentwise concatenates of Aijke (t) • An example of a member <l> n of the learning 
set is shown in Fig. 2b. 

The covariance T of the learning set was only implicitly formed as in (6). Since the matrix 
<l = ($1, . . . , is of size 1440 x 750, T would be of size 1440 x 1440. The calculation of the whole 
eigendecomposition of T would be a very burdensome task, especially when we only need a few 
principal eigenvectors. Fortunately, we can use the so-called orthogonal iteration method [12] to 
calculate the principal subspace in such a way that we do not even actually have to compute T. This 
method was then used to calculate 20 principal eigenvectors and the corresponding eigenvalues. 
The 20 largest eigenvalues of T are shown in Fig. 3a. The sections 7^(£), m = 1, 3, 5, 9 of the basis 
functions corresponding to the coefficient a 2 (t) are shown in Fig. 3b. 

Note that by construction of the learning set the times and the rates of change of the transitions 
as well as the initial and final states are assumed to be independent. 

3.2 Simulations 

To find a feasible dimension M of the basis we generated realizations with members from the 
learning set and formed the modified TVARLS estimates for different M. The choice M = 5 



6 



50 1 00 1 50 200 T 2T 3T 4T 5T 6T 

a) b) 

Figure 2: a) The extreme cases of sigmoids used in the construction 
of the learning set. b) An example of a member <£ n of the learning set 
showing the structure of the concatenated parameter evolutions. 




a) b) 

Figure 3: a) The 20 first eigenvalues of the covariance matrix T of the 
learning set <1> b) The sections of eigenvectors 7^(£), m = 1, 3, 5, 9 of the 
covariance matrix T corresponding to the parameter a 2 (t) . 
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Figure 4: a) The parameter processes a 2 (t) of three examples of 
the learning set (solid) and the corresponding estimates (dashed) based 
on single realizations, b) A true parameter process a 2 (£) (solid), the 
mean estimate (dashed) and standard deviation interval (dotted) of 20 
realizations. 



seemed to be a reasonable compromise between approximability and stability. Too low an order 
M means that the evolutions can not be adequately approximated by the basis while too high an 
order introduces instability to the estimates in sense that the realization-wise evolution estimates 
have non-acceptably high variance. 

As we are more interested in the difference between the true and estimated evolution of the 
prediction coefficients (mean deviation) than residual error variance we will focus our attention 
on the mean deviation of the parameter a 2 (t), the mean change of which is the greatest for the 
learning set. 

Three different evolutions of a 2 (£) and the corresponding estimates a 2 (£) from single realizations 
are shown in Fig. 4a. The initial and final states (and so also a 2 (l) and a 2 (T)) are equal in the 
three cases. As we can see, the estimates a 2 (t) can track the changes in a 2 (t) but tend to be biased 
towards the "center of the learning set" . In addition, the estimated rates of change are too low. 
These are features that can be explained by the nature of the eigenvectors as shown by Fig. 3a. The 
linear combinations of the 2 first eigenvectors have all approximately the same instant and slope 
of the transition. The features in 7fc(£) that allow for changes in the instant of the transition occur 
from eigenvectors 73 (t) on but all these have the same slope. The corresponding slope features 
occur in the eigenvectors from 79 (t) on but the features for estimating the instants of transition are 
still present. Thus the approximability of instants of transition of the evolutions would necessitate 
that we use M > 5. 

To test the variance of the estimates we calculated the means and variances of parameter 
evolutions in several cases of the learning set. A typical case is shown in Fig. 4b with true and 
mean evolutions and standard deviation intervals for 20 realizations. 

3.3 Tracking of ERS 

To evaluate the applicability of the method to real data we estimated the predictor evolutions for 
the data from which the initial and final predictor estimates were calculated, see Fig. la. 

Although we do not propose this method to be used in estimating the instant of an abrupt 
change in EEG, we conduct an experiment in which we try to evaluate the tracking capability of 
the estimates. Since there are now no true parameter processes, we can not compare the estimates 
to these. We take a sample x t of the measured EEG and form the estimates for x t+ T tv where T tr 
was allowed to extend somewhat beyond the modeling capability of the learning. 
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Figure 5: a) The ad hoc estimate of the change instand as a function 
of the shift of a single sample of EEC The slope is approximately linear 
(T tr ~ 0.75T tr + T 0 ) on the interval where the basis allows for changes in 
the parameters, b) An example of a modified evolutionary (amplitude) 
spectrum estimate based on the TVARLS model with var(e^) = 1. 



We estimate ad hoc the instant of change by the time T tr that is nearest to the mean of the 
estimated initial and final states, that is, a 2 (T tr ) & (22 (1) + a 2 (T))/2. These ad hoc estimates T tr 
as a function of the shifts T tr are shown in Fig. 5a. The relationship between these is approximately 
T tr ~ 0.75 • Tt r + To when T tr is compatible with the learning set. The reason that the slope is 
smaller than unity is that the employed set of basis has a clear bias to estimate the change instant 
towards the average of the change points in the learning set as explained in Section 3.2. This 
tendency is diminished when the order M is increased. We can also see that the degree of linearity 
diminishes when x t+ T tv is not compatible with the learning set. In Fig. 5b we show a typical 
evolutionary spectrum estimate (spectrogram) that is based on the time-frozen parameters of the 
estimated TVAR model. 

4 Discussion 

We have discussed a modification of the TVARLS scheme in which the individual parameter 
evolutions are not separately adjustable, but in which the mutual correlation of the parameters is 
taken into account. We have also presented a principal component type method for estimating such 
a basis for the parameter evolution that can optimally approximate the expectable evolutions. We 
have then applied the proposed method to the modeling of the ERD/ERS phenomenon of EEC 

The crucial problem of the method is the existence of a low dimensional subspace in which the 
expectable evolutions can be approximated with adequately small errors. If this is not the case, 
e.g. when the eigenvalues of T decay too slowly, this method should not be applied. The second 
problem is the formation of the learning set which may not always be as straightforward as in this 
case. 

It must be noted that the approximability of "marginal" members of the learning set with 
a basis is poorer than the approximability of more medial members. For this reason it could 
be recommended that the learning set should contain members that are slightly off the expected 
region. Otherwise it could be necessary to use a larger dimension M than with this trick if we 
wish to have a uniform approximability on the whole set. On the other hand, too large a learning 
set will yield a set of eigenvalues of T that will decay slowly, which means that to maintain the 
approximability we should again increase M which in turn decreases stability of the estimates. 

There are two previous approaches to estimate the ERD/ERS evolution. The first is a short 
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time Fourier transform -type (STFT) scheme in which discrete Fourier transforms are calculated 
from sliding windowed estimates of sample correlation [18]. In this method the window has to 
be made relatively short to maintain adequate time resolution and the variances of the time- 
frequency bins tend to become so large that the usefulness of the estimates suffers. We would thus 
be compelled to use averaging of the individual estimates. 

The other method is to pass the EEG realization- wise through bandpass filters, square the 
outputs and average over realizations [23]. Thus we obtain e.g. an estimate for the time-varying 
variance (power) of the alpha band. However, the bandpass filtering exhibits some drawbacks. 
The first is that the decay time of the filter impulse response (from 1 to 0.707) is approximately 
equal to the inverse of the width of the filter spectrum. If the band is 8-12 Hz, the decay time 
is about 0.25 s. In practice this means that the evolution estimates are effectively convolved with 
a window of this duration. The other problem in the filtering approach is that we are not able 
to track changes in the center frequency of a band within the passband. This implies that, e.g. 
if the alpha band "peak" moves towards the boundaries of the passband, but does not change in 
variance, the variance of the output of the filter decreases giving us a false indication. To overcome 
this problem we can divide the alpha band into sub-bands, but this will further decrease the time 
resolution. 

The proposed method does not necessitate ensemble averaging and it can thus track changes 
in single realizations and thus also trends through the ensemble. The frequency resolution is 
that of the AR models in general, but with the restrictions imposed by the subspace constraint. 
Some special information, such as the time- varying variances of the EEG bands can be directly 
accessed using the approximate spectral decomposition of the model, see [29] for the factorization 
of stationary parametric models. 

One further advantage of the proposed approach is that in some cases different kinds of 
deviations between parameter evolutions are portrayed in different basis vectors. It is thus possible 
that it is easier to draw inference from the coefficients of these basis vectors than from coefficients 
of generic basis. 
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